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Abstract 

This paper investigates the effect of drug treatment on the standard within-host HIV model, assuming 
that therapy occurs periodically. It is shown that eradication is possible under these periodic regimes, and 
we quantitatively characterize successful drugs or drug combinations, both theoretically and numerically. 
We also consider certain optimization problems, motivated for instance, by the fact that eradication 
should be achieved at acceptable toxicity levels to the patient. It turns out that these optimization 
problems can be simplified considerably, and this makes calculations of the optima a fairly straightforward 
task. All our results will be illustrated by means of numerical examples based on up-to-date knowledge 
of parameter values in the model. 

1 Introduction 

For the past two decades, within-host virus models describing the infection of HIV have played an important 
role in the understanding of this infamous retrovirus, and the ways in which it escapes not only the immune 
system, but also the various drugs that have been developed to suppress viral replication. Testing specific 
hypotheses based on clinical data is difficult since detection techniques of the virus are still far from accurate. 
This justifies the central role played by mathematical models in this area of research. 

An example of a question that has received considerable attention was whether drug treatment fails 
because of the pre-existence of drug-resistant strains, or by the emergence of resistant strains after initiation 
of drug therapy [4]. According to |14J . the former scenario is more likely. Nevertheless, the ability of the 
virus to mutate quickly, into forms which may be less sensitive to drugs has been, and continues to be, 
the focus of much attention, see recent contributions such as [2, 5 that study the behavior of multi-strain 
models. 

Other research has gravitated around the fact that the periodic regimen in which drugs are taken daily 
(or more frequently), puts a very high strain on the patient, calling for therapies minimizing the treatment 
burden [9], and also leading to investigations of the use of STI's (Structured Treatment Interruptions) 

[3 EES EES]- 

This paper revisits a by now classical model, often referred to as the standard model |13l 111] , which is a 
three-dimensional nonlinear ODE whose state consists of the concentrations of healthy CD4+ T cells (the 
targets of the HIV), infected T cells, and viral particles. Upon infection of a healthy T cell, one of the first 
orders of business is to make a copy of the viral RNA, using the enzyme reverse transcriptase. This step, 
which is error-prone and leads to mutations, can be blocked by a class of drugs called reverse transcriptase 
(RT) inhibitors. Once the viral copy has been produced, double stranded viral DNA integrates in the cell's 
nucleus as provirus. The usual gene expression now does the rest, and viral proteins are produced according 
to the genetic information encoded in the provirus. These proteins are assembled, mature and ultimately new 
viruses buds off from the infected cell's surface which go on to infect other T cells. During the maturation 
stage the protease enzyme is used to cleave long protein chains, and the so-called protease (P) inhibitors, 
are drugs that target this step. If effective, they give rise to defective virus. 

The purpose of this paper is to assess theoretically and quantitatively what the impact is of periodic 
drug treatment on the dynamic behavior of the standard model, and in particular to determine what it 
takes to get rid of the infection. Mathematically, we obtain a nonlinear periodic ODE, for which in general 
it is difficult to prove global stability and this explains why much research has traditionally resorted to 
simulations. Surprisingly though, solutions to the standard model ultimately arc bounded by solutions of 
a monotone system, as pointed out by d'Onofrio in [5], and this allows to conclude global stability for the 
nonlinear periodic model. 
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We will first consider a simple case, where only RT inhibitors are administered, and where it is assumed 
that the drug is of the bang-bang type, i.e. during a period of the treatment cycle, the drug is is either 
active or inactive. The drug is thus characterized by two parameters: its efficiency level when active, and 
the duration of the activity. A major role in our analysis is played by the spectral radius of a non-negative 
matrix (the fundamental matrix solution, evaluated over one period, of the linearization at the infection-free 
equilibrium), which is shown to possess expected monotonicity properties in terms of the two parameters 
that characterize the drug. Specifically, this spectral radius -which also controls the speed of convergence to 
the infection-free equilibrium- is lower when the drug is more potent or when it is active longer. Equivalently, 
convergence to the infection-free equilibrium is faster with a more potent drug, or a drug whose activity lasts 
longer. We will see that these results can be generalized to the case of P inhibitors, or to a mix of both RT 
and P inhibitors. This latter scenario reflects more closely the standard practice of administering cocktails 
of drugs to HIV infected patients. 

In reality, the efficiency of a drug is not of the bang-bang type. In fact, current research is investigating 
the effect of including pharmacokinetcs into the picture, and has revealed that the efficiency is a periodic 
signal with an initial steep rise right after drug intake, followed by a slower decay over a period, see the work 
of j7j [15] for detailed models. Therefore, we turn to this more general case, by approximating the efficiency 
by a more general piecewise constant periodic signal. It turns out that the previous results remain valid. 

Finally we turn to optimization problems that involve either maximizing the speed of convergence to the 
infection-free equilibrium while making sure that acceptable toxicity levels are not exceeded, or by minimizing 
toxicity levels, while making sure the speed of convergence does not fall below a certain threshold. 

All our results will be illustrated by means of numerical examples of within-host models whose parameters 
are chosen in accordance with current prevailing knowledge based on clinical data and extensive experimental 
evidence. Our results have the potential to suggest which drug, or which combination of drugs, are optimal 
for a given patient. They can also be used to explore the consequences of changing the treatment frequency. 
The investigation of the impact of periodic treatment cycles on multi-strain models, or the effect of STI's is 
the subject of ongoing research. 

Notation: For matrices A and B, < A, < A means that A is a (entry-wise) non-negative, positive 
matrix respectively, and A < B means that < B — A. A matrix is called quasi-positive if all its off-diagonal 
entries are non-negative. The spectral radius of a matrix A is defined as the largest modulus of all eigenvalues 
of A and will be denoted by p(A). 

2 Within-host HIV model with treatment 

We briefly recall the well-known standard model [T31 E] • Let 

T = f(T) - WT 
f* = kVT - (5T* 

V = Nf3T* - 7 V, (1) 

where T, T* , V denote the concentrations of healthy and infected T-cells, and virus particles respectively. 
All parameters are assumed to be positive. The parameters j3 and 7 are the death rates of infected T-cells 
and virus particles respectively. The infection is represented by a mass action term kVT, and N is the 
average number of virus particles budding off an infected T-ccll during its lifetime. The (net) growth rate 
of the uninfected T-cell population is given by the smooth function f(T) : R+ — > E, which is assumed to 
satisfy the following: 

3 T > : /(T)(T - T ) < for T ± T , and f'(T ) < 0. (2) 

We have chosen to make the class of allowable /(T)'s as large as possible, since the growth rate is hard 
to determine. In addition, most mathematical results apparently remain valid for this large class. Finally, 
we notice that the two most popular choices for /(T), namely a — bT for some positive a and 6, see 
and s + rT(l — T/T max ) for some positive s,r and T max , see |13j (here s is a source term modeling T cell 
production in the thymus and r and T max are the maximal per capita growth rate and carrying capacity 
respectively describing logistic growth of T cells), satisfy the preceding conditions. 
Since continuity of / implies that /(To) = 0, it is easy to see that 

T = (T , 0,0), 
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is an equilibrium of (JT|), and we will refer to it as the infection-free equilibrium. 

A second, positive equilibrium (corresponding to an infection) may exist if the following quantities are 
positive: 



Note that this is the case iff /(^) > 0, or equivalently by (|2j) that T — < T . In terms of the basic 
reproduction number 

_ kN 
Rq :— T , 
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existence of a positive equilibrium is therefore equivalent with 

1 < Ro, (4) 

which will be a standing assumption throughout the rest of this paper. Indeed, if we would assume that 
Ro < 1, it is known from [5] that the infection- free equilibrium Eq is globally asymptotically stable (GAS), 
and hence in this case the infection would always be cleared without treatment. 

We denote the positive equilibrium that corresponds to an infection by E = (T, T* , V). Linearization at 
Eq shows that it is unstable, and conditions on f(T) are known that guarantee that E is GAS (excluding 
of course initial conditions corresponding to a healthy, uninfected individual; these coincide with the T-axis, 
which is the stable manifold of E ). However, it is also possible that the model exhibits sustained oscillatory 
solutions which can be asymptotically stable. Regardless of the dynamical complexity of the solutions of the 
model, in general, if left untreated, the infection will persist within a patient. All these results follow from 
®. 

Obviously, the purpose of treatment is to clear the infection, hopefully by making Eq GAS by suitable 
modifications of model ([!} which reflect the effect of drugs. For the moment, we will only consider the effect 
of RT inhibitors, but P inhibitors will be included later. Using monotherapy based on RT inhibitors, model 
([I]) is modified to: 

f = f{T) — k(l — e{t))VT 
f * = k(l - e(t))VT - (3T* 

V = N/3T* - jV, (5) 

where e(t) £ [0, 1] is the (time-varying) drug efficiency of the RT inhibitors. The drug is not effective when 
e(t) = and 100% effective when e(t) — 1. Notice that Eq is still an equilibrium of the modified model 
regardless of the drug efficiency. 

Assuming that the efficiency is constant over time, we set e(t) = e £ (0, 1]. Then, to clear the infection, 
it suffices to choose e such that the modified basic reproduction number i?o(e) is less than 1, where 

«,(«) := «^T„. 
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Indeed, the results mentioned previously are applicable to this modified model, and they imply that if 
Ro(e) < 1, then E is GAS for Equivalently, if the efficiency e satisfies 

then treatment will be successful in this case. If the drug would be effective 100% so that e = 1, then 
treatment would always be successful. Current RT inhibitors clearly do not fit this profile. Moreover, in 
practice, the drug efficiency is not constant through time, and the main purpose of this paper is to investigate 
the quantitative consequences of this fact. 

3 Periodic drug efficiency 

We now make the assumption that e(i) is periodic: 

e(t) = e(t + t), for all t, 
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Figure 1: Periodic drug efficiency e(t). 



for some period r > 0. This is closer to reality where patients ideally adhere to a strict periodic treatment 
schedule, taking medication daily (r = 1 day) or twice a day (r = 0.5 day) for instance. 

The shape of e(t) over one period is determined by the (here unmodeled) pharmacokinetics, although 
coupling of the standard model with detailed pharmacokinetics models has been the subject of recent re- 
search, see for instance the work of [Zl [15], where it was shown that at least qualitatively, the graph of 
the periodic function e(t) is roughly like the one depicted in Figure [T| It is characterized by a quick rise 
of the efficiency to a peak value right after drug intake, followed by a slower decay. This is significantly 
different from the case where the efficiency is constant, the situation we described in the previous section. 
In pharmacokinetics, the efficiency e(t) is traditionally defined as 

<t)~ V{t) 



K + y{ty 

for some positive constant K. Here, y(t) is an output of a linear compartmental system 

z = Mz(t)+u(t), y(t) = z n (t), 

where M is a stable compartmental matrix (i.e. M is an n x n quasi-positive matrix whose eigenvalues are 
in the open left half plane), and the state components describe the concentrations of the drug in various 
compartments of the model (gut, blood, etc). Typically, y(t) is the concentration of an activated form of 
the drug inside the target cells. The input u(t) describes the (ideally periodic) drug intake signal and it is 
often modeled as a pulse. It is not difficult to show that for periodic u(t), the output y(t) will converge to 
a periodic signal with the same period, and this justifies to assume that e(t) is also periodic with the same 
period. 

Assuming a periodic efficiency e(t), let us start by linearizing system (|5| at the equilibrium Eq: 



x = B(t)x, (6) 



where 




-k(l - e(i))T N 
B{t)= ( -f3 fc(l-e(t))T 

-7 

It is well-known that the stability properties of the origin of ^ (and generically the local stability properties 
of the equilibrium Eq for system ([5])) are determined by the Floquet multipliers of ([6]). The block-triangular 
structure of B(t) implies that these are 

efl To > and A 2 ,A 3 , 

where A2 and A3 are the Floquet multipliers of the planar r-periodic system: 

z) - (jf t t(i -_r r °) ft) m 

In particular, since /'(To) < by ([2]), it follows that the three Floquet multipliers of system ^ are contained 
in the interior of the unit disk of the complex plane -which in turn implies that Eq is locally asymptotically 
stable for system ([5])- if | A2 1 , | A3 1 < 1. In fact, by a beautiful argument due to d'Onofrio in [8], it turns out 
that the same conditions imply the much stronger result of global asymptotic stability of Eq for system ([5]). 
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Figure 2: Periodic drug efficiency e(t) of the bang-bang type. 



Proposition 1. Let the Floquet multipliers of system |6| be contained in the interior of the open unit 
disk of the complex plane. Then E is GAS for system hence the infection is cleared. 

This result shows how relevant and important it is to determine the Floquet multipliers of system ^ . 
Unfortunately, for general functions e(t), this is a notoriously difficult task. Therefore, we will consider the 
simpler case where e(t) is piecewise constant, bearing in mind that piecewise constant functions are often 
good approximations to continuous functions. We will start with an even simpler case where e(t) is of the 
bang-bang type. 



3.1 Periodic drug efficiency of the bang-bang type 

We make the following simplifying assumption regarding the shape of the graph of the r-periodic function 
e(i), which is illustrated in Figure [2] 

where p £ (0, r) is the time duration during which the drug is supposed to be active with efficiency e <E [0, 1]. 
During the remaining part of the treatment period the drug is assumed to be totally inefficient. Clearly, 
this is a very crude way of approximating the more realistic shape of e(t) depicted in Figure [T| but some key 
properties are to be learned from this case, and they carry over to more general cases that describes reality 
better, as we will discover later. 

There are two possible parameters which can be varied in ([8]), namely e and p, and the purpose of the 
rest of this subsection is to investigate their effect on the Floquet multipliers of system (|7| with ^ . These 
Floquet multipliers arc the eigenvalues of the following matrix 



A{e)--=rl fc(1 " e)T °V (10) 



$(e ) p)=e( r - p ) jl ( V A ( e ' (9) 

where 

A(c\ — / 

sNP -1 

Since both A(e) and A(0) are quasi-positive matrices their matrix exponentials are non-negative matrices [^] 
Thus, <!>(e,p) is a non-negative matrix and by the Perron- Frobenius Theorem [3j its spectral radius p (<i>(e,p)) 
is an eigenvalue of <&(e,p). Thus, the Floquet multipliers of system ^ with ^ are contained in the interior 
of the unit disk of the complex plane if and only if p($(e,p)) < 1. This guarantees that the infection is 
cleared (globally) by Proposition [T| 

The following proposition -whose proof is deferred to the Appendix- reveals that p(<I>(e,p)) has the 
expected monotonicity properties: it decreases with e (more efficient treatment) and with p (drug is effective 
longer) . 

Proposition 2. Let e, e' € [0,1] andp,p' € [0, t]. Then the map (e,p) — > p(<fr(e,p)) is continuous, 

e<e', p^0^p($(e',p)) <p($(e,p)), (11) 



1 Proof: Let A be quasi-positive. Then B = A + al is a non-negative matrix for all sufficiently large values of a, implying 
that e tB is a non-negative matrix for all t > 0. But since e tA = e~ at e tB , the same conclusion holds for e tA . 
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Figure 3: Spectral radius of p (<&(e,p)) as a function of efficiency e and treatment duration p. The horizontal 
surface corresponds to p = 1 . 



and 

p<p', e?0=>p($(e,tf))<p(Q{e,p)). (12) 

Moreover, 

p($(e,0)) = p($(G,r)) = p(e rA( - 0) ) > 1 for all e 6 [0, 1] and aZZp S [0,r] (no treatment). (13) 

and 

p($(l,r)) = p (c tA(1) ) = max{e _/3T ,e" 7T } < 1 (constant, 100% effective treatment). (14) 

Since p($(e,p)) (provided it is less than 1) is a measure of how fast solutions of ^ with ^ approach 
E (at least locally near E ), this result may be interpreted as follows: 

Let the treatment be periodic, of the bang-bang type, and capable of clearing the infection. If it is more 
efficient, or lasts longer, then the infection is cleared more quickly. 

We illustrate Proposition [2] in Figures [3] and [4] The parameters used are taken from [TS] , and they are as 
follows: f(T) = a-bT with a = 10 4 ml" 1 day -1 and b = 0.01 day -1 (which implies that T = 10 6 ml" 1 ), 
k = 2.4 x 10" 8 ml day" 1 , (3 = 1 day" 1 , N = 3000, 7 = 23 day" 1 . The period of the treatment r is 1 day. 

Remark 1. This result can be modified to the situation in which P inhibitors are used for treatment instead 
of RT inhibitors. Model ([5| is then replaced by 

t = f(T) - kVT 
f * = kVT - (5T* 

V = N(l- e(t))0T* - 7V, (15) 



and matrix A(e) in (10) by 

A(o\ — ( 

K N(l-e)P - 7y 
With this notation and still using ([9]), Proposition [2] remains valid. 

Remark 2. Similar results can be stated to describe the situation in which combination therapy is used. This 
is the more commonly found therapy method where patients take a cocktail of both RT and P inhibitors. 
Model (|5| should then be replaced by 

f = f(T) - k(l - e RT (t))VT 
f* = k{l - e RT {t))VT - [3T* 

V = N(l-e P {t))/3T* -yV, (17) 
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Figure 4: Contour plot of the spectral radius p (<I>(e,p)) as a function of efficiency e and treatment duration 
p: p (<i>(e,p)) > 1 in red region and < 1 in green region. 
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Figure 5: r-periodic drug efficiency e(i) (dashed-dotted line) and a piecewise constant approximation. 



where 



e RT (t) = J eflT ' ' € [ °' PRT] , e P (t) = r P ' 1 G [ °' PP] 
RTK > l0, te(p RT ,T) py> lo, te(p P ,r) 



(18) 



denote the piecewise constant efficiencies of the RT and P inhibitors respectively. Finally matrix A(e) in 
(10 1 is replaced by 



A( e RT, ep) = 



N(l-ep)0 
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(19) 



With these notations and assuming without loss of generality that prt < Pp (if not, simply swap subscripts 
RT and P in the expression below), the spectral radius of the following matrix 

$(e flT , ep,p RT ,p P ) = e (T- P p)A(0,0) e { P P-p RT )A(0,e P ) e p RT A(.e RT ,e P ) 

is the key quantity. As expected, the spectral radius is increasing in each of its arguments epr, £p,Prt,Pp- 
We omit the proofs of these results as they are straightforward modifications of the proof of Proposition [2] 



3.2 General piecewise constant periodic drug efficiencies 

As mentioned earlier, in practice, the graph of the drug efficiency is not as shown in Figure [2] but rather 
as the dashed-dotted line in Figure [5] which can be approximated by a piecewise constant and r-periodic 
efficiency with several constant drug level efficiencies e% > €2 > ■ ■ ■ > e m > e m+ i during the respective 
intervals [p ,pi), [pi,Pa), ■ • ■ , bm-iiPm), bm -p m +i), where p Q := and p m+x := r for some m > 1. 
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Define (E,P) := (ei, . . .,e m+ x,px, . . . ,p m ) and let 

<f>(E,P) = e (P™-+i-P™) A ( e ™-+^ . . . e (Pi-Po)-A(ei) _ ^20) 

Similarly to Proposition [2] we find that the spectral radius of &(E,P) is increasing in each of its argu- 
ments. 

Proposition 3. Let m > 1, = po < p\ < P2 < ■ ■ ■ < p m < Pm+i = T and 1 = eo > e% > e2 > ■ ■ ■ > e m > 

P-m+i > e m+ 2 = 0. Then the map (E,P) — > p(<f>(E,P)) is continuous. In addition, 

ie{l,...,m + l} ande t+l <e i <e!' i <e i -R^p(^(E i ,P^ <p{$(E,P)), (21) 

and 

j G {1, . . . , m) and p 3 -_ x < ^ < pj < p i+1 p (* (tf, P./)) < p (* (£, P)) . (22) 

where Ei = (e^ e 2 , . . . , e^i, e-, e i+i , . . . , e m+1 ) and Pj = (pi,P2, ■ ■ ■ ,Pj-i,PpPj+i, ■ ■ ■ ,Pm)- 
The proof is deferred to the Appendix. 

4 Optimization problems 

In this section we return to the case of periodic efficiencies of the bang-bang type. What follows can easily 
be generalized to the case of more general, piecewise constant periodic efficiencies. As mentioned earlier, the 
purpose of treatment is to eradicate the infection by making Eq GAS for ([5| with In practice however, 
one would like to achieve this while the burden to the patient is as low as possible. Obviously, there are 
various ways to measure this burden. Let us list a couple of particular problems, assuming a r-periodic 
treatment schedule: 

1. Minimize p (<&(e,p)) subject to (e,p) £ [0, 1] x [0, r] and /J" e(t)dt — ep < c, for some fixed c € (0, r). 

2. Minimize JT e(t)dt = ep, subject to (e,p) £ [0,1] x [0, r] and p($(e,p)) < S, for some fixed 5 £ 
(p($(l,T)),l). 

In the first problem the spectral radius of <&(e,p) is minimized. As we mentioned before this spectral radius 
controls the rate of convergence to Eo (provided it is less than 1): the smaller the spectral radius, the faster 
solutions converge. In addition to minimizing the spectral radius, the burden on the patient should not 
exceed a specified upper bound c. Here, the burden to the patient is measured as the area under the graph 
of the efficiency e(t) over one period. The second problem on the other hand, concerns minimization of the 
patient's burden, subject to the condition that the spectral radius is less than a given bound 8 (assumed to 
be less than 1 so that convergence to E is guaranteed). 

Both problems fit in the larger classes of problems which we describe next. Let the maps F,G : [0, 1] x 
[0, t] — > M. be continuously differentiable with the following properties: 

F(0, 0) = G(0, 0) = 0, and VP, VG > 0, but ^ on [0, 1] x [0, r] \ {(0, 0)}. 

Now consider the more general optimization problems: 

Class I. Minimize p($(e,p)) subject to (e,p) £ [0, 1] x [0,r] and F(e,p) < c, for some fixed c > satisfying 

{(e,p)|P(e,p) = c}n[0,l] x [O,r]^0. 

Class II. Minimize G(e,p), subject to (e,p) £ [0, 1] x [0,r] and p($(e,p)) < 8, for some fixed 

*€(p(3(l,r)),l). 

The first two problems fit in this class for the choices F(e,p) = G(e,p) = ep. But it is clear that other 
choices could be of interest as well, for instance P(e,p) = ae qi + bp q2 , for some fixed q%, q2 > 1 and a, b > 0, 
or positive linear combinations of several of these functions. 

It turns out that both classes of optimization problems can be simplified thanks to Proposition [2] We will 
see shortly that the optimum appears on the boundary of the constraint set in both cases, which translates 

2 If i = l, then replace < by < in the right most inequality. If i = m + 1, replace < by < in the left most inequality. 
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Figure 6: Graph of p(&(e,p)) for points (e,p) with ep = 0.8 for e G [0.8,1]. Minimum for p is 0.665 and it 
is achieved for (e,p) — (1,0.8). 



into saying that the optimum occurs only if the patient's burden is the maximally allowed one (for problems 
in the first class), or that the spectral radius takes the largest allowed value (for problems in the second 
class) implying that convergence to Eo will be as slow as allowed. 

To be more precise, we claim that Class I and II optimization problems arc equivalent to Class III and 
IV problems respectively which are defined as follows: 

Class III. Minimize p ($(e,p)) subject to (e,p) G [0, 1] x [0, r] and F(e,p) — c, for some fixed c > 
satisfying that {(e,p)|F(e,p) = c} n [0, 1] x [0, r] ^ 0. 

Class IV. Minimize G(e,p), subject to (e,p) G [0, 1] x [0, r] and p(4>(e,p)) = 6, for some fixed 

*€(p(3(l,r)),l). 

Notice that the difference between Class I and III, and Class II and IV is in the constraint only (by replacing 
the inequality by an equality). In other words, the optimum of Class I and II problems occurs on the 
boundary of the constraint set. We show this equivalence for Class I and III problems. The argument to 
show equivalence of Class II and IV problems is very similar and omitted. Suppose that (e*,p*) is such that 
p(3>(e*,p*)) is minimal, while F(e* ,p*) < c. Notice that (e*,p*) ^ (0,0) since c > and F takes small 
positive values near (0, 0) in the rectangular region R := [0, 1] x [0, r] and p is strictly lower in those points. 
Also (e*,p*) ^ (1,t) since otherwise the level set {(e,p)\F(e,p) — c} does not intersect R, contrary to our 
assumption. 

If (e*,p*) is in the interior of R, then the point (e',p') :— (e*,p*) + sVF(e*,p*) is still in the interior of R 
with F(e',p') < c for small enough positive s, yet p (<&(e',p')) < p (<&(e*,p*)) by Proposition [2] contradicting 
minimality. The same argument applies if (e*,p*) = (0,p*) for some p* G (0,r) or if (e*,p*) = (e*,0) for 
some e* G (0, 1), since a perturbation of such a point in the direction of VF(e*,p*), results in a point which 
is still in R. If (e*,p*) = (l,p*) for some p* G (0, r) or if (e*,p*) — (e*,r) for some e* G (0,1), then a 
perturbation in the direction of V F(e* ,p*) could potentially result in a point outside R. To prevent this 
we perturb as follows for the case where (e*,p*) — (l,p*) (the argument when (e*,p*) — (e*,r) is similar 
and omitted): Let (e',p') — (l,p*) + (0, s). Then for sufficiently small and positive s, (e',p') is still on 
the boundary of R and F(e',p') < c, yet p(G>(e',p')) < p($(e,p)) by Proposition [2] a contradiction to 
minimality. 

5 Numerical examples 

Here we provide some examples of the optimization problems we just discussed. The model parameters used 
throughout this section are the ones chosen in subsection 3.1. 

Let us first minimize p ($(e,p)), see Figures [3] and [4] The constraint is that the burden to the patient, ep 
should not exceed 0.8. The minimum is 0.665 (which fortunately implies that with this treatment schedule 
the infection can be cleared successfully) and it is achieved at (e,p) = (1,0.8). In other words, the drug 
should be 100% efficient while it is active. This is illustrated in Figure [6] which depicts the spectral radius 
p($(e,0.8/e)). 

Let us see what happens when we modify the measure of the patient's burden to e 2 +p 2 , and demand 
that it should not exceed 1.2 2 . This time the minimal spectral radius is 0.866 (again implying that this 
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Figure 7: Graph of p(Q(e,p)) for points (e,p) with e 2 + p 2 = 1.2 2 for e <E [0.65, 1]. Minimum for p is 0.866 
and it is achieved for (e,p) = (0.868,0.829). 



p 




Figure 8: Level curves ep — c in red (c = 0.7,0.75,0.8,0.85,0.9, c increasing in NE direction). Level curve 
p($(e,p)) = 0.7 in blue. 

therapy will clear the infection) and it is achieved at (e,p) = (0.868,0.829). This is illustrated in Figure 

t which depicts the spectral radius p (<&(£, \/l-2 2 — e 2 )) . A striking difference between this schedule and 
e previous one, is that now the minimum is achieved in the interior of the rectangular parameter space 
[0, 1] x [0, 1], while previously it was achieved on the boundary. When the drug is active, it should therefore 
not be 100% efficient as before. 

Let us now consider minimization problems in which the patient's burden is minimized subject to a con- 
straint on the spectral radius, or equivalcntly, on the speed of convergence to the infection-free equilibrium. 
If the patient's burden is measured by ep, and if the spectral radius should not exceed 0.7, we find that the 
minimum is 0.782 and it occurs at (e,p) — (1, 0.782) which is on the boundary of [0, 1] x [0, 1] and requires 
that the drug is 100% effective when it is active. This is illustrated in Figure [8j where we depict a few level 
curves of ep, and the maximally allowable spectral radius p = 0.7. 

If we modify the measure of the burden to e 2 + p 2 (and still assuming the constraint that the spectral 
radius should not exceed 0.7), then the minimum is 1.493 and it occurs at (e,p) = (0.884,0.844) which is 
in the interior of the rectangular region [0, 1] x [0, 1]. This is illustrated in Figure [9] where we depict a few 
level curves e 2 + p 2 , and the maximally allowable spectral radius p — 0.7. 
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Figure 9: Level curves e 2 + p 2 = c in red (c = 1.4, 1.5, 1.6, 1.7, 1.8, c increasing in NE direction). Level curve 
p($(e,p)) = 0.7 in blue. 



Appendix 

Proof of Proposition [2] 

This proof hinges on the following standard facts: 

1. If A is quasi-positive and A < B but A ^ B, then < e tA < e tB but e tA ^ e* s for all i > 0. 

To see this, let a > be such that C = A + al > 0. Then setting D = B + al, we have that < C < D 
but C^B. Then e tC < e tD , but e tc ^ e*- for all i > 0. It follows that e tA < e tB but e tA ^ e tB for 
all i > 0. 

2. For all t > 0, e tA ( ) > if e ^ 1, while c tA ^ > (but not > 0) if e = 1. 

3. HA > and B > has no zero row or zero column, then AB > and BA > 0. This is true in 
particular when B = e tc for t > and C a quasi-positive matrix because of Fact 1 and the fact that 
matrix exponentials are invertible. 

4. If < A < B but B ^ A, then p{A) < p(B), see Corollary 1.5 in Chapter 2 in [3]. 

Continuity of the map (e,p) — > p(Q(e,p)) follows from the definition ^ of p and the fact that the spectral 
radius of any matrix is continuous in terms of its entries. 
Let < e < e' < 1 and p ^ 0. Then: 

A(e') < A(e) and A(e') ^ A(e) 
=> < e p ^ e ') < e pA(e) and e pA(e,) ^ e pA ^ by Facts 1 and 2 

< e^-P^ e^ e '> < e ( T -f^(°) e^ e ) and e^^ e? A ( e '> f e^ e > by Facts 1 and 3 

and invertibility of matrix exponentials 
=> < <J>(e',p) < <I>(e,p) and $(e',p) ^ $(e,p) 
=> p($(e' : p)) < p($(e,p)) by Fact 4. 



This result remains valid if e' = 1 because /?(<I>(e,p)) is continuous. This establishes (111 
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Let < p < p' < t and e^O. Then: 



A(e) < A(0) and A(e) ^ A(0) 

< e (p'- P )A(e) < e (p'~p)A(0) and e (p'-p)A(e) ^ e (p'-p)A(0) by Fact j 
< e (p'-P) A (e) e pA(e) < e (p'-p)A(0) e pA(e) &n( j 

e ( P '-p)A(e) e pA(e) ^ e ( P '- P )A(o) e pA(e) by Fact j and invertibility of exponentials 

< q(t-p')M0) e (P~P)A(e) e pA(e) < e (-r-p')A(0) g(p'-p)A(0) e pA(e) and 
e (T-p')A(0) e (p'~p)A(e) e pA(e) ^ e (r-p')A(0) e (p'-p)A(0) e pA(e) by 3 and 

invertibility of exponentials 

< <&(e,p') < $(e,p) and <&(e,p') ^ $(e,p) 

p($(e,p')) < P(He,p)) by Fact 4. 



This remains valid if p' — t because p(<&{e,p)) is continuous. This establishes (12l. 

Finally, it follows from our standing assumption Q, that the determinant of A(0) is negative. Thus ^4(0) 



has a positive eigenvalue which implies (13 1. Also, (14 1 is immediate from (10) 



Proof of Proposition [3] 

The same facts as in the proof of Proposition [2] will be used. 

Continuity of the map (E,P) — > p(Q(E,P)) follows from the definition ^ of p and the fact that the 
spectral radius of any matrix is continuous in terms of its entries. 

Fix i G {1, . . . ,m + 1} and let a < e- < 1. Then: 

A(e'A < A( ei ) and A(e[) + Afa) 

< eto-**- 1 )^ < e (Pi - pi ~ l)Aiei) and e {p ^- l)A{ < ] ^ e (p *- p4 - l)A(ei) by Facts 1 and 2 
=S> < $ (Ei, p) < $ (£■, P) and $ (.E;, P~) ^ $ (P, P) by Facts 1 and 3 

and invertibility of matrix exponentials 
=> p(<f>(Ei,p)) <p ($(P, P)) by Fact 4. 



This result remains valid if i = 1 and = 1 because p(^(E 1 P)) is continuous. This establishes (21 1. 
Fix j € {1, . . . , m} and let < pj < p'j < r. Since e^+i < ej, we have that 

A(e 3 ) < A{e j+1 ) and A( ej ) ^ A{e j+1 ) 

=S> < e ( ^^ )j4(e ^ < e (p i- p > )A(e ' +l) and e (p '^ p ^ A{e ^ ± e ( ^" Pj)A(ej+l) by Fact 1 
=>■ < e^^'^* 6 ^ e( pj '~ Pj_1 )' 4 ( e; '' < e (Pj-P3)^( e 3+i) etPj-Pj-O^t 13 ^ an( j 

e (^-p 3 )A(e 3 ) e ( Pj - Pj -i)A( e3 ) _^ e (p;- P3 )A(e 3+1 ) gfe-ft-iMfe) by Fact : and invertibi i ity of exponentials 

=4> < e ( - Pj + 1 ~ P '^ A( - ej + 1 ^ e^Pj^P^-^M^i) < e (Pj + i-Pj)A(e j + 1 ) e (p j -p j -i)A(e j ) &nd 

e (P 3 + l-^M(e 3 + l) e (p;-p 3 -l)A(e 3 ) ^ gCpj+i-pjJXCe^Ogfe-Pi-OA^) by Facts X and 3 and 

invertibility of exponentials 
=> < $(P, Pj) < $(P, P) and $(P, P,-) ^ $(P, P) by Facts 1 and 3 and invertibility of exponentials 

=> p(<J>(P,P,)) <p($(P,P)) by Fact 4. 



This establishes (22) 
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